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^. 

^T^^ In recent decades a lot of research has been done on the numerical solution of the time- 

dependent Schrodinger equation. On the one hand, some of the proposed numerical 
methods do not need any kind of matrix inversion, but source terms cannot be easily 
C^i ' implemented into this schemes; on the other, some methods involving matrix inver- 

\^ sion can implement source terms in a natural way, but are not easy to implement into 

\J , some computational software programs widely used by non-experts in programming (e.g. 

xl i Mathematica) . We present a simple method to solve the time-dependent Schrodinger 

Si-J ' equation by using a standard Crank-Nicholson method together with a Cayley's form for 

(_^ ' the finite-difference representation of evolution operator. Here, such standard numerical 

j/~v , scheme has been simplified by inverting analytically the matrix of the evolution operator 

^~v ■ in position representation. The analytical inversion of the N X N matrix let us easily 

^~v ' and fully implement the numerical method, with or without source terms, into Mathe- 

matica or even into any numerical computing language or computational software used 
for scientific computing. 

- Keywords: Schrodinger equation; Finite-difference methods; Numerical simulation; 

J^ ' Mathematica 6.0. 

d ' PACS Nos.: 01.50.H-, 02.70.Bf, 02.10.Ud, 02.60.Cb, 03.65.Ge 

1. Introduction 

One of the main arguments used to explain why quantum mechanics is not easily 
accessible to most of the students attending for first time a quantum mechanics 
course, is that the physical situations of quantum mechanics are not everyday life 
phenomena and that can be approached only through abstract mathematics. The 
significant differences between the classical and quantum physics, makes difficult to 
our macroscopically-trained minds to imagine what is happening in a physical situ- 
ation at the quantum regime. On the other hand, usually in the classical mechanics 
courses as a general problem-solving strategy is suggested first to draw a sket ch or 
diagram that represents the physics of the problem under consideration.^^ Such 
strategy is contradicted in the case of quantum mechanical problems where in al- 
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most any introductory quantum mechanics t extb ook can be found statements as: 
do not try to imagine the physical situation!^'^ Told to avoid visualization, the 
students fall into misunderstandings because the lack of a mental picture leads to 
inefficient problem solving'- . 

Fror n th e previous discussion and as has been considered for some other 
authors}^'-^ it is necessary to implement visualization techniques that can improve 
both, understanding of, and problem-solving in quantum mechanics. In the present 
paper we concentrate particularly in particle propagation methods. Along this line, 

Q 

since the pioneering work by Feit et al.° some efforts has been done aiming for a 
comprehensive and easy implementation (even for beginners in the field) of numeri- 
cal methods to solve the Schrodinger equation (see for instance Ref.lHland references 
therein). Yet, it should be noted that nowadays there are many programs that are 
able to do such simulations) ^^ ^ but in most cases, those are numerical codes 
that require more than a basic knowledge of programming to be implemented. 

Writing a simple program that can be fully implemented by the student al- 
lows a better understanding of the phenomenon, opens the possibility of treating 
a wider range of problems and allows the student to have a pleasant first con- 
tact with programming. In this paper we describe a numerical integration method 
for the time-dependent Schrodinger equation and the implementation of absorb- 
ing bou ndary c onditions and source terms into this scheme. Unlike most existing 
methods) I I this method is easy to implement, versatile, very accurate, and can be 
implemented without including any kind of matrix inversion package, which allows 
a fully implementation on a mathematical package such as Mathematica, providing 
easily visualizable results of the evolving system. 

The plan of this paper is as follows. In section [5] we present a short description 
of the technique and its numerical implementation. We follow by presenting the 
formula for the inverse of a non-symmetrical tridiagonal Jacobian matrix and we 
show how to introduce the inverse matrix formulas into the numerical scheme. In 
section [3] we continue relating it to the absorbing boundary conditions to finally 
apply it for the case of the source term in section |4l We end this paper by compar- 
ing our numerical results with the analytical solution for the Schrodinger equation 
with source term to finally compare the analytical and numerical results for the 
transmission probability of a finite potential barrier. 



2. The Method 

Let us consider the Schrodinger equation in atomic units, i.e. m — h — 1, 

d 

with the Hamiltonian given by 

H{x,t)^~^-^ + V{x,t). (2) 
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The idea is to compute the time evolution of the wave function ip{x,t) for t > 
to, given an initial state ?/'(x,to)- We start by dividing the time interval into n 
subintervals of equal length At = {t ~ to)/n, and use an implicit Crank- Nicholson 
integrator scheme L!^ to propagate the wave function from one time step to the next 
one. 

The formal solution to Eq. ([T]) could be expressed in terms of the time evolution 
operator as, 

i^{x,t)^e~'"'il;{x,0). (3) 

The effective time evolution operator U for one discrete time step Ai, can be ex- 
pressed using Cayley's form for the finite-difference representation of e^*^*, which 
is a combination of a fully implicit and a fully explicit method, 



l-i^H{x,t) 
l + i^H{x,t)' 



U{t + M,t)= ^ , ^1,J:^: . (4) 



Such representation oilA is second-order accurate in space and time and also unitary. 
The integration scheme for the wave function then reads 

1 + '-^H{x, t) j ij{x, i + At) = M - ^^(^' *)) '^(^' *)• (^) 

The wave function can be expanded on a discrete lattice as 

N 

V'(a:,i„)=^V>, (6) 

i=i 

where "0" = '4'i^j: tn) is the value of the wave function at the position Xj of the jth 
lattice site at time in = ^o + nAt, with a grid basis 

^ ^ r 1, Xj - i Ax <x<Xj + i Ax; 
'' \ 0, otherwise. 

Here Ax = (xmax — a^min)/-^, with Xmax and Xmin the boundaries of the finite grid. 
Using the finite-difference representation for the kinetic part of the 
hamiltonian,L!^ we have 

1 ± fn) .,.,.« . ,.; . f (- ■^?---^;^-- . vr^ (8) 

with V" = V{xj,tn). By introducing tp^ — (?/)", ...,?/;",..., ip^), the lattice represen- 
tation of Eq. ([5]) finally reads 

i/;"+i = D-^DiV^"" , (9) 

where we define 

Di=(l-^i/)=(l-S), D2=[l + ^i/)=(l + S), (10) 



December 2, 2010 1:11 WSPC/INSTRUCTION FILE ws-ijmpc 



4 F. L. Dubeibe 



with S = ^^H. The matrix product can be rewritten as 

Dj^^Di = (1 + S)-i(l - S) = 2T>^^ - 1, 



(11) 



then, the wave packet evolution is achieved just by inverting the matrix D2. 

For the case of time independent potentials, the explicit N x N representation 
of Di and D2 reads 



Di = 



/71 


a 
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a 


72 
a 


a 

73 


a 
a 


IN- 


1 a 




\ 








a 


In 


J 



D, 



/a -a 
-a ^2 —a 
-a £,3 -a 



\ 



-a £n^i -a 
-a £n J 



with 



a 



iAt 
4:Ax' 



iAt ( 1 



r, 7j = 1 - ft-, e^ = 1 + ft, and /3, = ^- 



Aa;2 



y, 



(12) 



(13) 



2.1. Inverse of a Tridiagonal Matrix 

Let us consider the N x N nonsingular tridiagonal matrix D 

I ai bi \ 

ci 02 62 
C2 as ^3 



D 



V 



(14) 



CN-2 flTV-l &Ar-l 

cn-1 aw y 

UsmaniLtSlgave an elegant and concise formula for the inverse of the tridiagonal 
matrijifj: 






(15) 



where 9i satisfy the recurrence relation 

6lj = aj6'i_i - 6i_iCi_i6'i_2, for i^2,...,N, (16) 

with initial conditions ^q = 1 and 61 = fli, and 0j satisfy the recurrence relation 

4)^ ^ ai(f)t+i - biCi(j)i+2, for i = iV - 1, . . . , 1, (17) 

with initial conditions 0jv+i = 1, <pN ~ ajv, and On = detD. 



*A few typos and misprints from the original paper were corrected. 
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Using the last procedure, we obtain a simplified formula for the inverse of the 
matrix D2, 

(D2X-.i = rf,, =(-l)'+-'(-a)IJ"-*l0,_i,^,+i/0A,, if 1<J. (18) 

Due to the tridiagonal symmetric nature of D2, the inverse satisfies, dji — dtj. The 
recurrence relations are given by 



and 

4>i = 

with 9q = I, 61 = ai. 
Finally, the elenw 

by 



where Sij is the Kronecker delta. 

In what follows we consider two standard systems belonging to the class of time 
independent potentials: The finite square potential well and the finite square barrier. 
In both cases Dirichlet boundary conditions are assumed. These boundary condi- 
tions may cause unwanted reflections, therefore, we have to perform the numerical 
calculations on a sufficiently large bounded interval, placing the impinging particle 
far away of the numerical boundary and restricting the time interval such that the 
reflections do not affect the solution in the region of interest. 

2.2. Finite Square Potential Well and Finite Square Barrier 

In order to test the method and to observe the effect of the Dirichlet boundary 
conditions, we consider a Gaussian wave packet 



— Ci^i-l ~ Ct 9i-2, 


for i = 2,...,iV, 


(19) 


^i4>i+i - a^0i+2, 


for i = iV- 1,...,1, 


(20) 


(/)Ar+i = 1, and 07V 


= aN- 




nts of the matrix product E = Ti2^Ti\ = 2D^^ - 


- 1, are given 


(E),, =2d,J-5,J 


for z, j = l,...,iV, 


(21) 



V'(a;,0) = //^-exp 

V ^0^^ 



IPq[X - Xq) ;^--2 — 



2oi 



(22) 



initially centered at xq, with average momentum po and initial width ctq, which 
moves into the region of a short range potential defined as 

Vix-] - [ ^^'/^' ~^b < ^ < ^b; ^23) 

^ ' \0, otherwise. ^ ' 

This is a potential barrier (plus signed) or a potential well (minus signed) whose 
height or depth, respectively, equals the average energy Po/2 of the Gaussian wave 
packet.L!^ In Fig. [1] and Fig. [2] we show the Gaussian wave packet scattering from 
the finite square barrier and the finite square well respectively. 

As can be seen in Figs. [T] and [21 there exists a strong back refiection effect due 
to the Dirichlet boundary conditions in both cases. The wave packet behaves like 
inside a large infinite square well of length L — 2a;,„ax- If the wave packet spreads 
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quickly, any reflected portion of the wave will then interfere with the portion of 
the incident wave, giving rise to a non-physical interference pattern. This situation 
imposes limitations on the choice of the input parameters, e.g. xq and uo in Eq. 
([22)) must be chosen so that ip{—L/2,0) and ip{L/2,0) are essentially zero at least 
at the beginning t — to. 




f 


t=3.6 


f\ Ll_ 


— - 



9 t=4.3 



iL 




Fig. 1. Gaussian wave-packet scattering from a finite square potential barrier. The initial condi- 
tions are xq = —10, itq = 1, Xb = 2, to = and po = 7. The left and right borders of the domain 
are Xniin = —20 and imax = 20 respectively. The size of every lattice in the grid is Ax = 0.04 with 
N = 1000 discrete lattices in all the spatial domain and the time step is At = 0.002. With the 
given parameters the matrix elements are calculated from Eq. I I13I1 . The parameter t = tg -|- n At 
denotes the time of each configuration. 



For the square barrier case. Fig. [TJ a fraction of the wave packet is captured 
by the barrier and remains trapped for a period which is longer than the time of 
transmission through the barrier. The captured piece of wave packet bounces back 
and forth between the barrier walls with a small amount of probability escaping in 
each collision till the entire packet escapes. With the present numerical scheme the 
dynamic evolution of the trapped wave-packet can be easily observed at each time 
step. 



3. Method with Absorbing Boundary Conditions (ABC) 

The numerical solutions of the time-dependent Schrodinger equation provide us in- 
sight into the dynamics of quantum mechanical systems. However, as was discussed 
in the previous section in practical calculations the area of computation must be 
limited to a finite grid because of the finite capacity of the computer memories. This 
finite grid produces undesirable reflections at the artificial boundaries of the area of 
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Fig. 2. Gaussian wave-packet scattering from a finite square potential well. The same parameters 
as in Fig. (T] 



computation. To minimize this artificial effect we implement in the present section 
the so called absorbing boundary conditions. These are local boundary conditions 
that approximate the one way wave equation of a wave function. 

Let us set the following equation as the starting point of the discussion 



., 9 , , , 



-|l^ + n.))^(.,i). 



(24) 



In orde r to obtain the formulas for the absorbing boundary conditions, following 
Shibata}^ we consider the special solutions i^{x,t) = exp{—i{Lot — kx)), these are 
states of definite energy E satisfying the dispersion relation, 



hk = ±y^2m{huj-V). 



(25) 



The absorbing boundary conditions must be designed to satisfy the dispersion 
relation given by the plus signed Eq. (|25|) at the boundary Xmax and the minus 
signed at the boundary a^min- However, function (I25p is not rational and cannot 
be converted into a partial differential equation, nonetheless, this relation can be 
linearly approximated by 



hk = gi{huj -V)+g2, 



(26) 



with 



91 =± 



^/2ma2 — \/2mai 
a2 — ai 



92 =± 



a2\/2mai — OL\\J2ma2 



a2 — Oil 



(27) 
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The correspondence of d/dt <^ —ioj and d/dx <^ ik leads us to rewrite Eq. (^7)) 
into the partial differential equation 

d , . / . I 



ifi^i'ix^t) 



-ih- 



gidx 



V - 



51/ 



(28) 



Now we outhne how to incorporate the ABC into the lattice representation of 
the wave function (with h — m = l). The idea is to replace the differential equation 
for the boundary components f/'Ar and ?/)" of the state vector -0" . As was discussed 
by Paul et aLj^in order to obtain an accurate expression for the derivative at the 
borders of the grid is convenient to introduce an intermediate point x between the 
last two points of each side of the grid, then for example, at the right hand side the 
wave function must be replaced by 

V'(x,i)~-[i/;(a;w,i)+V(xjv-i,i)]. (29) 

In the grid representation, the finite-difference equation for the right and left sides 
reads 



2At 



—Ml^' 



V'^-i 



rN-rN-i) 



(jiAa; 



irN - 1 



N-l) 



+ l{^-f)^'^N + rN-i), 



(30) 



and 



2Ai 



(vr^+vr^ 



i^^ - V^n 



giAx 



(V'2 - ^l) 



f^)ir. 



-r,. 



(31) 



respectively. The equations ([30|) and (|3T|) allow a straightforward incorporation into 
the matrix representation. 

The new matrices T)i2 are given by 

a j3 a —a ^3 —a 



Di 



with 



a 7Ar-i a 

V3 m) 



,D2 



-a ^N~i -a 

m m J 



(32) 



-'^^-2At' 






i z 1 


/ 

V - 


92 


~ 2At giAx ' 2 


51 


i i 1 




.92 


~ 2At giAx ' 2 


51 



(33) 
(34) 

(35) 
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The main cause of artificial back reflection for plane waves in the presence of 
the above boundary conditions comes from the approximate nature of the finite 
difference evaluation. Clearly, these approximations become better decreasing the 
grid spacing Ax. 

In the present case, the matrix elements in (1141) are given by ai — ajv = ?72, 
bi = cn-1 = 111, bj = —a for j = 2, ..., A^ — 1 and Cj — —a for j — 1, ..., N — 2. The 
elements in (D2)ij^ = dij are given by 

dii = (t)2/dN, 

di, = {-l)^+^fl2{-ay'-^^<t>j+,/0M, for J = 2,...,N 

d.j = (-l)'+^(-«)l-'"*le.-i0,+i/eAr, for i,j^2,...,N with i<j 

dij = {-iy+^i-ay'-^\9j^i(j)i+i/9N, for i = 2,...,N-l, j = l,...,N with i>j 

dAT, = (-l)^+^77i(-a)l^-J"-i|0,_i/0^, for 3^l,...,N^l. 



with 



70 



= 1, 



^1 = m- 

e^^^^e,^i-a^e,-2, for z-3,...,A^-1 



and 



(/•Af+l = 1, 
4>N = V2, 

01 =^i0i+i - a^(?!)i+2, for i = A^-2,...,2 
01 = »y202 +a?7i03- 

On the other hand, the non-symmetric character of the matrices Di^2 does not let 
us to write the matrix product E = Ti^^Tii as simple as in Eq. (|TT|) . The new 
components of the product are 

E.a = r]idii + adi2, for i^l,...,N 

Ei2 = mdii +I2di2 + adt3, for i = l,...,iV 

Eij ^ adij-i +^jdij + adij+i for i = 1, ..., N,j ^ 3, ...,N ~ 2 

EiN-i = adiN-2 + iN-idiN-i + mdiN for i = 1, ..., N 

E,N = ad^N-1 + VidiN for i = 1, ..., N 



It should be noted that the parameters in Eq. (fT3|) are the same, but now these 
parameters are defined between 2 and iV — 1. 
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3.1. Finite Square Potential Well and Finite Square Barrier with 
ABC 

In order to observe the effect of the ABC at the boundaries of the grid, we use the 
same examples as before: The square well and the square barrier with the Gaussian 
wave packet Eq. (P^ as test particle. The potential function is defined as in Eq. 
([23l) . In Fig. [3] and Fig. IH we show a wave packet scattering off a square well and a 
square barrier respectively. 



a 


t=0.0 
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t=0.7 




c 


A 


t=1.4 









d t=2.2 




<= t=2.9 




f t-3.6 








Fig. 3. Gaussian wave-packet scattering from a finite square potential well in the presence of 
ABC. In this case we choose cti = 24 and a2 = 25, the other initial conditions and parameters are 
the same as in Fig. [2] 



As can be seen from the Figs. [3] and SI the boundary conditions effectively 
reduce the non-physical reflections of the impinging wave packet at the boundary 
of the computation area. The clear presentation of the resonances let us study its 
dynamical evolution in a detailed way. 



4. Method with Absorbing Boundary Conditions and Source Term 

As a last example illustrating the validity and effectiveness of our simplified method, 
we consider the presence of a source term. The equation of motion now reads 



d 

i-—ip{x, t) — Hix, t)^p(x, t) -\- Sit) exp{—iujt)5ix) 
ot 



(36) 



where S{t) = 5*0 [1 — cxp(— t/AT)]. This function provides a smooth evolution of 
the source term towards the desired final value S'(i -> oo) — 5o. In order to have an 
analytical result, we consider the stationary solutions of Eq. (pS)) for the particular 
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Fig. 4. Gaussian wave-packet scattering from a finite square potential barrier in the presence of 
ABC. The same parameters as in Fig.[Tl 



case V{x) = with 5(t) = 5o. Introducing the Fourier transformed wave function 
il^{q,t) = j ex.'p{iqx)ip(x,t)dq and considering the previous conditions, the Eq. (155)) 
takes the form 



dt 2m 

This equation admits solutions of the form 

2^0 



tp{q,t) — 5*0 exp(— iwi). 



(37) 



i'iq,t) 



fc2 



exp(— iwt) 



where k^ = 2lu. Transforming back to the configuration space, the solution is given 

by 



(38) 



ip(x,t) = — ;- exp(i/c|a;|) exp(— iwi), 
ik 

showing that the source emits in bo th d irection a monochromatic wave Eq. 

As it was shown by Paul et. al.}^^ working with a grid representation of the 

wave function, it is convenient to approximate the S function by 



R(x) == -r^[e(.T + Aa;/2) - e(x - Aa:/2)], 
Ax 



(39) 



where Q is the Heaviside step function. With this approximation, the error scales 
quadratically with the size of the grid Aa; and becomes negligible for reasonable 
small values of Ax. The implementation of the source term at position Xji in the 
grid representation reads as 
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where 5jj' = I ii j = j' and otherwise. In the presence of the source term, Eq. 
dH) is given by 



where the components of 6" are defined as 



(40) 



(41) 



The numerical implementation is the same as in the previous section with the only 
difference that we have to construct the new vector b and subtracts it as was indi- 
cated in Eq. (J40| . i.e. we can use the results given above for the analytical inverse 
D^ and for the product E = Dg^ Di. 



4.1. Plane Waves with Constant Amplitude 

Here we consider the case V{x) — and S{t) = Sq in which the exact solution was 



given as 



c 

ip{x,t) — — - exp(jfc|x|) exp(— iwi). 
ik 



(42) 



In Fig. [5] we compare the exact Eq. (|42|) and the numerical results for this case. 
The agreement between the numerical and exact result suggest that the method 
is sufficiently accurate and stable. It should be noted the excellent behavior of the 
Absorbing Boundary Conditions, which, regardless that the source is filling the 
numerical region, the numerical evolution simulates an open domain. 
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Fig. 5. Refi/i] as a function of distance x, for the analytical (dashed line) and numerical (dots) 
solutions to the Schrodinger equation with source term. The initial conditions are xq = 0, to = 0: 
ai = 12 a2 = 13, So = 5, lu = (po~92)/9i and po = 5- The left and right boundaries of the domain 



-10 and Xn 



10 respectively. The size of every lattice in the grid is Ax = 0.02 



with N = 1000 discrete lattices in all the spatial domain and the time step is At = 0.001. With 
the given parameters the matrix elements for l|32|l are calculated from Eq. l |13|l . The given times 
are for the variable t = to -|- nAf . 
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4.2. Plane Waves with time dependent Amplitude 

The case of plane waves with time-dependent amphtude is of physical interest, 
because the idea of an initially empty waveguide that is gradually filled with matter 
waves corresponds to the experimental realization of a reservoir located at x — xq. 
For propagation times t ^ Ai, the calculation converges toward a flat density that 
corresponds to the stationary plane waves at the source amplitude S = Sq. The time 
evolution of the probability density during the increase of the source amplitude is 
displayed in Fig. [6l 
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Fig. 6. |i/i| as a function of distance z, for the numerical solution to the equation with time- 
dependent source term. 



The transmission through a potential barrier has been a model of great relevance 
from a pedagogical point of view, as discussed in almost every textbook on quan- 
tum mechanics. The exact solution to this problem is usually obtained by assuming 
that a plane wave impinges on the barrier from the left. However, the comparison 
of the analytical result for the transmission and/or reflection coefficients with the 
numerical one is not an easy task, because the usual calculation in terms of gaus- 
sian wave-packets gives us an average of the analytically calculated transmission 
coefficients. This fact can be understood taking in to account that a Gaussian wave 
packet can be viewed as a superposition of plane waves with different momentum. 
Then, the correct way to compare the transmission and/or reflection coefficients in 
this case is to solve the numerical problem with a source term emitting plane-waves. 
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In Fig. [71 we compare the analytical transmission coefficients for the case of a finite 
potential barrier with the predicted by our numerical method with time-dependent 
amplitude. Also in this case we find an excellent agreement with the exact results. 



0.6 



0.2 - 




Analytical 
Numerical 



10 

Vn 



Fig. 7. Comparison of the analytical and numerical results for the transmission probability T of 
a finite potential barrier, as a function of the potential height Vb. The same parameters as in Fig. 

m 



5. Summary 

We have described a numerical integration method for the time-dependent 
Schrodinger equation with an without source terms. In particular, we consider the 
case of scattering systems, in which the Dirichlet boundary conditions produces un- 
desired refiections. To solve this problem we have introduced absorbing boundary 
conditions (ABC) into this integration scheme. The numerical integration was done 
using the Crank-Nicholson method together with a Cayley's form for the finite- 
difference representation of evolution operator which produces an stable, unitary, 
and second-order accurate in space and time method. On replacing the Hamiltonian 
by its finite-difference approximation, the problem reduces to a complex tridiago- 
nal system. We have simplified the numerical scheme by inverting analytically the 
matrices by means of the Usmani's formula for Jacobian matrices. The analytical 
inversion of the matrices decrease the computational effort and let us fully and easily 
implement the method into Mathematica or even into any scientific computational 
software. This numerical method can be used for arbitrary potential shapes even in 
the presence of the source term, which does not limit its use to elementary applica- 
tions required in teaching quantum mechanics. The formalism discussed here may 
be extended in a straightforward way to time-dependent potentials, in which the 
matrices vary in every time step. Finally, we have compared the results of our mod- 
ified numerical method with the analytical solution for the transmission probability 
of a finite potential barrier, where we find an excellent level of agreement. 



December 2, 2010 1:11 WSPC/INSTRUCTION FILE ws-ijmpc 



Solving the time-dependent Schrodinger equation in Mathematica 6.0 15 

Acknowledgments 

We enjoyed fruitful discussions with K. Rapedius, M. Hartung, L. A. Pachon, T. 
Dittrich, K. Richter, and C. Viviescas. Financial support from Volkswagen Foun- 
dation (grant 1/78235), Universidad Nacional de Colombia in the program Becas 
para Estudiantes Sobresalientes de Posgrado, Colciencias, and the ALECOL pro- 
gram of the German Academic Exchange Service DAAD is gratefully acknowledged. 
We thank for the hospitality extended to us by the MPI for the Physics of Com- 
plex Systems, Dresden, University of Technology, Kaiserslautern, and University of 
Regensburg, where part of this work has been carried out. 

References 

1. P. J. Nolan, Fundamentals of college physics (Wm. C. Brown Publishers, 1995). 

2. J. W. Jewett, R. A. Serway, Physics for scientists and engineers with modern physics, 
Seventh edition (Cengage Learning EMEA, 2007). 

3. D. Halliday, R. Resnick, J. Walker, Fundamentals of physics, Volume 2, Fourth edition 
(Wiley, 1993) page 1173. 

4. D. F. Styer, The strange world of quantum mechanics (Cambridge University Press, 
2000) page 115. 

5. B. S. Ambrose, P. S. Shaffer, R. N. Steinberg, and L. C. McDermott, Am. J. Phys. 
67, 146 (1999). 

6. N. S. Rebello, D. A. ZoUman, Visual quantum mechanics: progress report (Kansas 
State University, 1997) 

7. B. Thaller, Advanced Visual Quantum Mechanics (Springer, India, Pvt. Ltd., 2008); 
B. Thaller, Visual Quantum Mechanics (Springer, India, Pvt. Ltd., 2009s) 

8. M. D. Feit, J. A. Fleck, Jr., and A. Steiger, J. Comput. Phys. 47, 412 (1982). 

9. M. H. Degani and M. Z. Maialle, J. Comput. Theor. Nanosci. 7, 454 (2010) 

10. R. M. M. Mattheij, S. W. Rienstra, J. H. M. ten Thije Boonkkamp, Partial differen- 
tial equations: modeling, analysis, computation, Volume 10 of SI AM monographs on 
mathematical modeling and computation (SIAM, 2005). 

11. J. C. Strikwerda, Finite difference schemes and partial differential equations (SIAM, 
2004). 

12. K. Liu and A. Wagner, The chemical dynamics and kinetics of small radicals, Volume 
2 (World Scientific, 1995). 

13. R. Kosloff, Ann. Rev. Phys. Chem. 45, 145 (1994). 

14. T. Fabcic, J. Main and G. Wunner, J. Chem. Phys. 128, 044116 (2008). 

15. W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flanery Numerical Recipes 
in C++ (Cambridge University Press, Cambridge, 2002). 

16. W. A. Ames, Numerical Methods for Partial Differential Equations (Academic, New 
York, 1977). 

17. A. Goldberg, H. M. Schey and J. L. Schwartz, Am. J. Phys. 35, 177 (1967) 

18. R. Usmani, Linear Algebra Appl, 212 413 (1994). 

19. S. M. Blinder, Am. J. Phys. 36(6) 525 (1968). 

20. T. Shibata, Phys. Rev. B. 43 6760 (1991). 

21. T. Paul, M. Hartung, K. Richter and P. Schlagheck, Phys. Rev. A., 76, 063605 (2007). 



